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I .  INTRODUCTION 


One  duty  of  the  United  States  Coast  Guard  is  the 
interdiction  of  drug  smugglers.  When  Coast  Guard  personnel 
board  a  suspect  vessel  to  search  for  contraband,  a  decision 
must  be  made  regarding  the  amount  of  time  allocated  for  the 
investigation.  If  a  careful,  but  lengthy,  search  is 
conducted,  there  is  a  high  probability  of  finding 


contraband  if  it  is  present, 
investigation  incurs  a  penalty. 


However ,  a  lengthy 


While  Coast  Guard 


personnel  are  searching  one  vessel,  other  shipping,  which 
could  have  contraband  onboard,  is  passing  through  the 
patrol  area  without  being  examined.  At  the  other  extreme, 
if  the  suspect  vessel  is  examined  in  a  cursory  manner,  the 
Coast  Guard  Vessel  (CGV)  will  be  able  to  stop  and  search 
more  shipping  but  there  is  a  higher  probability  of  not 
discovering  the  contraband,  if  present,  because  of  the 
short  investigation  time.  The  purpose  of  this  thesis  is  to 
develop  a  renewal-reward  model  to  determine  the  optimum 
investigation  time. 


II.  RENEWAL-REWARD  MODEL 


Consider  the  following  scenario.  All  shipping 

encountered  by  a  patrolling  CGV  is  stopped  and  searched  to 
determine  if  contraband  is  onboard.  After  leaving  port,  a 
CGV  searcher  a  time  Si  until  a  vessel  is  sighted.  A 
detected  vessel  has  a  probability  Pb  of  being  Bad  (having 
contraband  onboard)  and  PB  ■  1  -  Pb  of  being  Good  (having 
no  contraband  onboard) .  A  Bad  vessel  contains  the  random 
number  J  =  j(l,2,3,...)  units  of  illicit  substance  with 

probability  bj .  The  time  required  to  find  the  first  unit 
of  the  contraband,  given  that  J  units  are  present,  is 
modelled  here  as  the  minimum  of  J  independent  identically 
distributed  times  having  distribution  function  Fz  (t) ;  this 
model  is  illustrative  only,  and  may  be  altered  in  various 
realistic  directions.  If  time  T  is  required  to  locate  an 
incriminating  unit,  then 

P(T>T|J-j|  =  [l-Fz(t)]J  ;  j-1,2,3,... 

so,  upon  removing  the  condition, 

•• 

1-Fr  (t)  -  P(T>t)  »  S  [1-Fz  (t )  ]-*  bj  ;  j-1,2,3,...  (2.1) 

H 

Adopt  the  following  decision  rule:  Establish  a 

predetermined  investigation  time  L.  When  a  vessel  is 
detected,  it  is  stopped  and  searched.  If  no  contraband  is 


discovered  in  the  investigation  time  L,  the  vessel  is 
released  and  the  CGV  resumes  patrol.  If  any  contraband  is 
discovered  before  the  end  of  the  investigation  time  L,  the 
vessel  is  detained  and  escorted  to  port  for  further  search 
and  investigation.  A  time  period  D  is  required  to  escort 
the  vessel  to  base.  The  patrol  cycle  ends  when  the 
escorted  vessel  arrives  at  the  base. 

The  long  run  average  reward  per  cycle,  R,  can  be 
calculated  in  the  following  manner  [Ross:  pp.  279-294]. 

E[Rc] 

R  -  -  (2.2) 

E[C] 

where:  E[Rc]  is  the  expected  reward  for  a  patrol 

cycle . 

E[C]  is  the  expected  duration  of  a  patrol 
cycle . 

The  reward  associated  with  an  apprehension  can  be  defined 
in  various  ways.  Three  different  reward  criteria  will  be 
examined.  For  each  of  the  three  cases,  we  will  determine 
the  optimal  investigation  time  L  which  maximizes  the  long 
run  average  cycle  reward. 

A.  EXPECTED  DURATION  OF  PATROL  CYCLE:  E[C] 

He  must  distinguish  between  the  following  two  cases  for 
the  investigation  time  I : 

la  -  an  investigation  time  that  results  in  the  release  of 
the  vessel. 

Ib  -  an  investigation  time  that  results  in  the  detention 
of  the  vessel. 


r 
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Where 


la  =  L  with  probability  PB  +  Pb Ft  (L)  =  a 

lb  *  T  with  probability  Pb  Ft  (L)  *  o’ 


The  length  of  a  CGV  patrol  cycle  can  be  represented  as: 


C  * 


Si  +  lb  +  D  With  probability  a(L) 
Si  +  la  +  C'  With  probability  a(L) 


where  C'  has  the  same  unconditional  distribution  as  C,  for 
if  the  first  investigation  results  in  a  release  of  the 
vessel  being  searched,  the  process  re-starts  (regenerates) . 

The  expected  length  of  a  patrol  cycle  can  be  expressed 
as : 


or 


E [C]  *  E[S]  +  (E [lb ]  +  E [D] ) a  +  (L  +  EtC*])a 

E f Si  +  E [Pi a  +  E fib  1  a  ±  La 


E[C] 


which  can  be  written  as 


(2.3) 


(2.4) 


1  -  a 


E  [C]  -  E[D]  +  E[Ib]  +  +  L| 

a  oi 


Solving  for  the  term  E[Ibl  yields 

r  tFi  <dt> 

EIIbl  -  -  F,  [L) 
Integration  by  parts  shows  that 


(2.5) 


(2.6) 


Ft  ( t )  dt  -  Ft  (L)  L  +  tF(dt) 


i: 


(2.7) 


f 


r, 

-i 

3 

i\ 


'A 

1 


t  jM 

V 

4 


Pjj 


s 

s 

s 

V 


3 

■faMacaaMMiMaai^^ 


4 


Hence 


3 


a 


E[Ib] 


G  Ft 


(t)dt  -  LFt  (L) 


Ft  (L) 


(2.8) 


Therefore 


EEC] 


E[S]  +  E  [D]  Pb  Ft  (L)  +  Pb  Fr(t)dt  +  LP9 


Pb  Ft  <L) 


(2.9) 


which  reduces  to 


EEC]  =  E[D]  +  EES]  +  j0  Ft  (t)dt  +  p„  L 

Pb  Ft  (L)  Ft  (L)  Pb  Ft  ( L ) 


(2.10) 


B.  REWARD  CRITERIA  1 


Suppose  the  CGV  is  rewarded  for  making  an  arrest. 


Since  only  one  arrest  is  made  per  cycle. 


E[Rc]  -  1 


(2.11) 


and  from  equation  2.2 


EEC] 


(2.12) 


as 


Consider  the  following  example 


Fz  ( t )  =  1  -  e'Ut 

,  P  =  0 

.08 

S(t)  =  1  -  e~Xlt 

,  Xi  = 

1.1 

D  ( t )  =  1  -  e'Xzt 

,  Xa  = 

0.3 

bj  =  (1  -  P)Pj_1 

,  P  =  0 

. 8 ;  j  =  1,2,... 

Pb  =0.2 

From 

equation  2.1 

-  M 

Ft  (t)  =  2E  (e-“‘  )i 

H 

(l-p)pj-i 

e-»*  (1-p) 

=  (2.13) 

1  -  pe-wt 

and 

Ft  (t)  =  1  -  ¥7(t)  = 

1  -  e-  m  t 

(2.14) 

1  - 

From 

equation  2.10,  the 

expected 

duration  of  the  patrol 

cycle  is 


and  from  equation  2.12,  the  long  run  average  reward  per 
cycle  can  be  calculated. 


1-e-uL 


l-fie- u 1 


Xi  X2  Pb 


(2.16) 


(  \  /  1-p  \  /  1-pe-ML  \ 

Pb  -  Xi  +  X*  +  Xi  X2  Pb  -  Ln  -  )  +  Pg  LXi  X2 

V  1-pe-uL/  V  &U  y  V  1~P  / 


C.  REWARD  CRITERIA  2 

The  reward  assigned  to  the  CGV  is  the  amount  of  drugs 
confiscated . 


E  [Rc  ]  =  E  [  J  |  T  S  L] 


(2.17) 


The  probability  that  a  vessel  being  searched  contains  J 
j  units  of  drugs  and  that  the  drugs  are  discovered  is 


P{ (T£L) ~(J-j) j 


(l-Fi(L)H  bj  (2.18) 


and  thus 


P{TsM«Ft(L)  =  2:  Tl  -  (1-Fz  (L)  ) J  lbj  (2.19) 

H  L  J 


jl  -  ( 1-Fz 


(L)  )J  bj 


Pf (J-j) I (TSL) I  * -  (2.20) 

Ft  (L) 


E[J|T£L] 


S  j  [l  -  (1-Fz  ( L )  H  "I  bj  (2.21) 

i=l  L  1 


using  the  parameters  from  the  Example  in  Model  1, 


2  j  [l  -  (e-«L  (l-p)pj-i 


E [ J | F£L]  =  - 

( 1-e"  )/  (1_ge-uL } 

which  can  be  evaluated  as  in  the  previous  example  as 


(2.22) 


E[J|T£L]  = 


(1-6)*)  -  (e'“  ] 
(i-e,[(V  u-»)  - 


(2.23) 


using  equation  2.2,  the  long  run  average  reward  is 


E [ J  TSL1 


E[C] 


[(—  )-  (— - )1 

L  \(1-P)V  \  (1-Pe-»L  )  V  J 

] 

L  \  (1-p)  /  \  (l-pe-««-  )  /  J 


(2.24) 


i-  (— 

Xi  \  Pu/  \  1-P  /  P«L 

N  /1~e~ML  \  p„  \ 

\  1-P e-“L/  \  1-Pe-«L  /  \1-Pe-ML/ 


l-e- ml 


D.  REWARD  CRITERIA  3 

The  reward  examined  in  model  three  is  the  difference 
between  the  quantity  of  drugs  confiscated  and  the  quantity 


of  drugs  that  escape  detection.  The  expected  reward  in  a 
cycle  is 

Rc  -  E [ J | TSL]  -  E [M] 

where  M  is  defined  as  the  quantity  of  drugs  that  escape 
detection  in  a  cycle. 

f 

Jib  +  Jd 
Ji  «  +  Jo  +  M 


where : 


Jib  «  Drug  quantity  passing  through  patrol  area  while  CGV 
searches  a  Bad  ship 

Jo  -  Drug  quantity  passing  through  patrol  area  while  CGV 
escorts  a  Bad  ship  to  shore 

Jib  *  Drug  quantity  passing  through  patrol  area  while  CGV 
searches  a  Good  ship 

Jo  *  Drug  quantity  onboard  a  ship  classified  as  Good 
that  is  actually  Bad 

If  conditional  expectations  are  taken  assuming  that  the  Bad 
vessels  pass  through  the  region  according  to  a  Poisson 
process  with  rate  Xi Pb ,  then 

E [M]  -  Xi  PbB[J]  (Btlb  ]+E[D]  )PoFt  (L)  +  (Xi  PbE[J]E[Ig  ]+E[M]  ) 

(Pa  +Pb  Ft*  (L)  )  +E  [Jo  ]  Pb  Ft"  (L)  (2.25) 


With  probability  PbFr(L) 

With  probability  PB  +  PbFt  (L) 


$1 
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*'i3 


sg 
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53 


>51 
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At  Pb  E  [  J] 


Ft  (L)dt-LFL  I  +E  [D J  Pb  Ft  (L)+L  Pb +Pb  Ft  (L 


PbFT  (L) 


E  t  Jo  ]  Pb  Ft  <L) 


PbFT  (L) 


(2.26 


where  E[Jo]  is  defined  as 


E [ J I T>L]  = 


E  j(l-Fz  (L)  Hbj 
j=l _ 

Ft  (L) 


Using  the  distributions  defined  in  the  Example  for  Model  1 


1-Be->*L 


( 1-B )  /  1-e~  \l 

\  1-Be->«L  / 


l-e-  u  l 


1-p  /  A2(l-Be->*L)  Pb 


(2.27) 


(l-B)e-»>- 

( l-pe~ <* L  )  (l-e-»«L) 


and  the  long  run  average  reward  can  be  expressed  as 


R  =  Equation  [2.211  -  Equation  [2.271  (2.28) 

Equation  [2.15] 


V.V.V.V 
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III.  RESULTS 


The  optimum  investigation  time  was  determined  for  the 
above  example  for  each  reward  criteria  using  two  methods. 
The  first  method  consisted  of  writing  a  Monte  Carlo 
simulation  of  the  patrol  cycle  using  the  distributions  of 
the  example  and  running  the  simulation  using  various  values 
of  the  investigation  time  L.  Statistics  were  gathered 
during  the  simulation  allowing  the  calculation  of  the  long 
run  average  reward  R.  The  investigation  time  was  varied 
from  zero  to  six  hours  in  increments  of  two  tenths  of  an 
hour.  The  optimum  investigation  time  was  determined  by 
graphing  the  long  run  average  reward  as  a  function  of 
investigation  time  and  finding  the  value  of  L  which 
maximized  R.  Each  cycle  was  replicated  20,000  times.  A 
detailed  discussion  of  the  simulation  can  be  found  in 
Appendix  A. 

THe  second  method  consisted  of  writing  a  computer 
program  for  the  three  equations  representing  the  long  run 
average  reward  for  the  three  reward  criteria  for  the 
example  and  solving  the  equations  for  various  values  of  the 
investigation  time  L.  As  in  the  case  of  the  simulation, 
the  value  of  L  was  varied  from  zero  to  six  hours  in 
increments  of  two  tenths  of  an  hour.  The  optimum 
investigation  time  was  again  determined  by  graphing  the 


long  run  average  reward  as  a  function  of  investigation  time 
and  finding  the  value  of  L  which  maximized  R. 

Figures  1,2,  and  3  contain  the  results  of  both  the 
simulation  and  numerical  solution  using  the  distributions 
of  the  examples  for  reward  criteria  1,2,  and  3 
respectively.  The  solid  line  represents  the  analytical 
solution  and  the  circles  are  the  simulation  results.  The 
maximum  long  run  average  reward  using  criteria  1  can  be 
achieved  by  using  an  investigation  time  between  1.4  and  2.2 
hours.  The  maximum  long  run  average  reward  using  criteria 
2  and  3  can  be  achieved  using  an  investigation  time  between 
1.2  and  1.4  hours.  It  is  interesting  to  note  that,  using 
the  distributions  presented  in  these  examples,  and 
investigation  time  exists  that  maximizes  all  three  reward 
criteria  simultaneously.  This  occurs  at  1.4  hours  for  the 
input  distributions.  We  do  not  know  that  this  state  of 
affairs  will  persist. 
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IV.  SENSITIVITY  ANALYSIS 


In  order  to  check  the  robustness  of  the  model  with 
respect  to  the  underlying  assumptions,  the  simulation  was 
run  exactly  as  before  with  the  exception  that  F2  (t)  is  now 
assumed  to  have  a  lognormal  distribution  instead  of  an 
exponential  distribution.  Each  cycle  was  replicated  20,000 
times.  Three  separate  cases  are  examined.  In  case  A,  the 
distribution  function  Fz  (t)  has  the  same  mean  and  variance 
as  the  exponential  distribution  used  in  the  previous 
examples.  In  case  B,  the  distribution  function  F2  (t)  has 
the  same  mean  but  twice  the  variance  as  the  exponential 
distribution  used  in  the  previous  examples.  Finally,  in 
case  C,  the  distribution  function  Fz  (t)  has  the  same  mean 
but  four  times  the  variance  as  the  exponential  distribution 
used  in  the  previous  examples. 

Figures  4,5,  and  6  contain  a  comparison  between  the 
exponential  base  case  and  the  lognormal  cases  mentioned 
above  for  reward  criteria  1,2,  and  3  respectively.  It  is 
readily  apparent  that  the  results  for  the  exponential  base 
case  and  lognormal  case  A,  with  the  same  mean  and  variance, 
are  different.  Furthermore,  within  the  lognormal  family  of 
curves,  it  can  be  seen  by  examining  case  B  and  C  that  the 
results  vary  significantly  as  the  variance  increases. 


LEGEND 

EXPONENTIAL  -  BASE  CASE 
o  LOGNORMAL  -  CASE  A 
a  LOGNORMAL  -  CASE  B 
+  LOGNORMAL  -  CASE  C 


LEGEND 

EXPONENTIAL  -  BASE  CASE 


LOGNORMAL  -  CASE  A 
LOGNORMAL  -  CASE  B 
LOGNORMAL  -  CASE  C 
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Figure  5:  Reward  Criteria  2  Sensitivity  Analysis 
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Figure  6:  Reward  Criteria  3  Sensitivity  Analysis 


In  order  to  better  understand  this  behavior,  it  is 
helpful  to  compare  the  quantiles  of  T,  the  time  required  to 
locate  an  incriminating  unit  given  that  J  units  are 
present,  for  the  four  distributions. 


iuantile 

Base  Case 

Case  A 

Case  B 

Case  C 

0.10 

0.242 

1.604 

0.908 

0.431 

0.25 

0.852 

2.502 

1.467 

1.681 

0.50 

2.465 

3.945 

2.856 

1.681 

0.75 

6.504 

6.858 

5.360 

3.909 

0.90 

13.488 

12.155 

10.333 

11.245 

It  can  be  seen  that  some  distributions,  in  particular 
the  lognormal,  cases  a,  b,  and  c,  have  a  far  greater 
likelihood  of  producing  large  values  of  T  relative  to  the 
others,  such  as  the  exponential  distribution  used  in  the 
examples.  Since  the  time  required  to  find  the  first  unit 
of  drugs  greatly  influences  the  cycle  length,  distributions 
generating  larger  values  of  T  will  produce  significantly 
different  results.  To  further  verify  the  accuracy  of  the 
distributions  produced  by  the  simulation  and  the  quantiles 
listed  above.  Appendix  B  contains  the  calculated  quantiles 
for  lognormal  case  A  produced  using  analytical  distribution 
theory. 
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V.  CONCLUSIONS 


It  has  been  demonstrated  that  a  renewal-reward  approach 
to  modelling  the  Coast  Guard  drug  interdiction  process  is 
feasible  and  that  it  is  possible  to  determine  an  optimal 
investigation  time.  By  considering  several  different 
reward  criteria  in  the  model,  it  is  possible  to  quantify 
the  costs  with  respect  to  the  rest  of  the  reward  criteria 
when  one  criteria  is  selected  and  used  to  arrive  at  an 
optimum  investigation  time. 

Of  the  two  methods  presented  to  obtain  the  optimum 
investigation  time,  the  simulation  approach  is  the  most 
flexible.  The  example  distributions  and  assumptions  used 
in  this  thesis  were  selected  because  they  allowed  an 
analytical  solution  to  be  found.  This  allowed  a  comparison 
of  results  between  the  simulation  and  numerical  solutions, 
thus  verifying  the  simulation.  As  the  model  assumptions 
are  changed  and  different  probability  distributions 
incorporated,  the  simulation  can  be  easily  modified  to 
reflect  these  changes  whereas  an  analytical  solution  may  no 
longer  be  possible. 

As  the  model  is  currently  formulated,  its  usefulness  is 
questionable  due  to  the  sensitivity  to  the  underlying 
distribution  of  the  time  required  to  find  drugs  onboard  a 
vessel  given  drugs  are  present.  As  demonstrated  in  the 
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sensitivity  analysis  even  within  a  given  family  of 
distributions,  the  model  is  sensitive  to  the  distribution 
parameters.  In  order  for  the  model  to  produce  realistic 
results,  rigorous  data  analysis  must  be  conducted  to 
properly  identify  this  distribution. 

The  basic  service  system  model  may  have  wider 
application,  e.g.,  to  combat  modelling,  where  it  may  be 
desirable  to  investigate  a  potential  target  to  estimate  its 
value  before  committing  limited  resources.  An  adaption  of 
the  model  may  also  be  of  help  in  allocating  resources  for 
mineral  exploration. 
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APPENDIX  A 


The  simulation  is  written  in  Fortran  77.  It  consists 
of  a  main  program  and  seven  subroutines.  Uniform  [0,1] 
random  numbers  are  provided  by  calling  procedure  GGUBFS 
located  in  the  IBM  IMSL  single  precision  library.  All  real 
variables  are  computed  using  double  precision  to  minimize 
rounding  error.  Output  is  directed  to  three  separate 
units.  Unit  2  contains  the  calculated  rewards  for  each 
increment  of  investigation  time  using  reward  criteria  1,  2, 
and  3.  Unit  3  contains  the  cycle  length,  drug  quantity 
confiscated,  and  quantity  of  drugs  missed  for  each 
increment  of  investigation  time.  Unit  4  contains  detailed 
information  regarding  the  cycle  for  each  increment  of 
investigation  time.  This  information  includes  the  total 
number  of  ships  searched,  the  number  of  ships  that  are 
good,  bad,  bad  but  declared  good,  and  bad  identified  as 
bad.  Also  included  in  unit  4  output  is  the  time  required 
to  discover  the  drugs  on  a  ship  declared  bad,  the  quantity 
of  drugs  onboard  ships  passing  through  the  area  while  the 
CGV  is  searching  a  ship,  the  quantity  of  drugs  missed  due 
to  short  investigation  time,  and  the  quantity  of  drugs 
confiscated  on  a  bad  ship. 

The  main  program  controls  the  starting  and  final 
investigation  time,  investigation  time  increment  width,  and 


simulated.  Cycle  averages  are  computed  and  output  directed 
to  the  three  units  discussed  above. 

Subroutine  CYCLE  simulates  one  patrol  cycle  and  records 
all  the  relevant  statistics  during  the  cycle.  Subroutine 
STIME  generates  the  random  search  times  required  to  find  a 
ship  from  an  exponential  distribution.  Subroutine  DTIME 
generates  the  random  times  required  to  escort  a  ship  back 
to  base  from  an  exponential  distribution.  Subroutine 
DRGQTY  generates  the  random  quantity  of  drugs  on  a  BAD  ship 
from  a  geometric  distribution.  Subroutine  CLASS  determines 
the  classification  of  a  vessel  based  on  the  deterministic 
value  Pi>  .  Subroutine  OPCOST  determines  the  quantity  of 
drugs  missed  onboard  other  ships  passing  through  the  area 
while  the  CGV  is  investigating  the  current  ship. 
Subroutine  RTIME  generates  the  random  times  required  to 
find  the  first  unit  of  drugs  given  that  j  units  are 
present. 

Four  versions  of  this  subroutine  are  listed;  one  for 
the  exponential  distribution  and  the  other  three  for  the 
lognormal  case  A,  B,  and  C  distributions. 
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FILEi  DPSIM 


FORTRAN  A1 


THIS  PROGRAM  SIMULATES  PATROL  CYCLES  FOR  VARIOUS  LENGTHS  OF 
INVESTIGATION  TIMES  AND  CALCULATES  THE  VALUES  FOR  THE  THREE 
REWARD  CRITERIA.  FUNCTION  GGUBFS  FROM  THE  IMSLSP  LIBRARY  IS 
USED  TO  GENERATE  UNIFORM  (0,1)  NUMBERS. 

INTEGER  I , K,MAXT >  NCYCLE, Z 

REALX8  L,C,J,M,INC,EC(500),EJ(500), EMC 500 ) , ITI ME (500), TO END, 

5  R1(500),R2(500),R3(500),TS,GS, BS, BSG, BSB, EIB, 0P,MISS, EIG, 

6  CATCH, ES, ED, BJ, DSEED 
DSEED  3  995317. DO 
MAXT  3  4 
NCYCLE  =  100 
INC  -  0.20 
L  3  0 

TOEND  =  MAXT/INC 
DO  200  K  3  1, TOEND 
L  3  L  +  INC 
C  3  0 
J  3  0 
M  3  0 

PRINT*, 'COMPUTING  L  3  ' , L 
DO  100  I  3  1, NCYCLE 

CALL  CYCLEC DSEED, L,C, J,M,TS,GS,8S, BSG, BSB, EIB, OP, MISS, 
B  CATCH, ES, ED, BJ, EIG) 

100  CONTINUE 

Z=NCYCL£ 

WRI TEC 4 , x )  L 

WRITEC 4^523)  TS/Z, GS/Z, BS/Z, BSG/Z, BSB/Z 
WRITE(4,524)  EIB/Z, BJ/Z, OP/Z,MISS/Z,CATCH/Z 
WRITE(4,525)  EIG/Z, ES/Z, ED/Z 
WRITEC4,*) 

523  FORMAT (1X.5F9. 4) 

524  FORMAT ( IX, 5F9 . 4) 

525  F0RMAT(1X,3F9 .4) 

TS  =0 
BS  =0 
BSG30 
BSB=0 
GS  30 
EIB=0 
MISS=0 
CATCH=0 
ES  =0 
ED  =0 
BJ  =0 
OP  =0 
EIG=0 

ITIMECK)  3  L 
EC(K)  3  C/NCYCLE 
EJCK)  3  J/NCYCLE 
EM(K)  3  M/NCYCLE 
Rl(K)  3  l/ECCK) 

R2(K)  3  EJ(K)/EC(K) 

R3(K)  3  (EJ(K)-EM(K) )/EC(K) 

WRITE(2,19)  ITIME(K) ,R1 (K),R2(K),R3(K) 

WRITE(3, 19 )  ITIMECK), ECCK) , EJ (K) , EM(K) 

19  FORMAT  ( IX, F6 . 2, 3F9 .3) 

200  CONTINUE 

STOP  'PROGRAM  COMPLETE' 

END 


SUBROUTINE  CYCLECDX, LSTAR, CSTAR, JSTAR, MSTAR, TS, GS, BS, BSG, BSB, 
a  EIB, OP, MISS, CATCH,  ES, ED,  BJ, EIG). 

REAL *8  LSTAR, CSTAR, JSTAR, MSTAR, LAMDA1 , LAMDA2, MU, BETA, S,R,QTY, 
a  D, LEAK, TIME, PBAD,TS,GS, BSG, BSB, EIB, OP, MISS, ES, ED, BJ,BS, CATCH, 
a  EIG, DX 
INTEGER  NC 
CHARACTERX4  TYPE 
LAMDA1  3  1.10 
LAMDA2  3  0.3 
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FORTRAN  A1 


C\ 


20 


21 


MU  =  0.08 
BETA  =0.8 
PBAD  =  0.2 
TS=TS+1 

CALL  STIME(DX,LAMDA1,S) 

ES=ES+S 

CALL  CLASSCDX, PBAD, TYPE) 

IF<  TYPE  .EQ.  'GOOD*)  THEN 
EIG  =  EIG+LSTAR 
GS=GS+1 
QTY=0 

TIME  =  LSTAR 

CSTAR  =  CSTAR  +  TIME  +  S 

CALL  OPCOSTCDX.LAMDAl, TIME,  BETA, PBAD, LEAK) 

MSTAR  =  MSTAR  +  LEAK 

0P=0P+LEAK 

GOTO  20 

ELSE  IF(  TYPE  .EQ.  'BAD')  THEN 
BS  =  BS+1 

CALL  DRGQTY( DX, BETA, QTY) 

BJ=BJ+QTY 

CALL  RTIME(DX,MU,QTY,R) 

IF(R  .GT.  LSTAR)  THEN 
EIG=EIG+LSTAR 
BSG=BSG+1 
TIME  =  LSTAR 

CSTAR  =  CSTAR  +  TIME  +  S 

CALL  OPCOST ( DX, LAMDA1 , TIME, BETA, PBAD, LEAK) 

MSTAR  =  MSTAR  +  LEAK  +  QTY 

OP=OP+LEAK 

MISS=MISS+QTY 

GOTO  20 

ELSEIFCR  .LE.  LSTAR)  THEN 
EIB=EIB+R 
BSB=BSB+1 

CALL  DTIMECDX, LAMDA2,  D) 

ED=ED+D 

TIME  =  R  +  D 

CSTAR  =  CSTAR  ♦  TIME  ♦  S 

JSTAR  =  JSTAR  +  QTY 

CALL  OPCOST(DX,LAMDAl, TIME, BETA, PBAD, LEAK) 
MSTAR  =  MSTAR  +  LEAK 
OP=OP+LEAK 
CATCH=CATCH+QTY 
GOTO  21 


ELSE 


STOP  'ERROR  1' 


ENDIF 
ELSE 

STOP  'ERROR 
ENDIF 
RETURN 
END 


SUBROUTINE  STIME< DRSEED, INPUT, OUTPUT) 
REALX8  INPUT, OUTPUT, RV,DRSEED 
RV=GGUBFS( DRSEED) 

OUTPUT  =  LOG(RV)/( -INPUT) 

RETURN 

END 


SUBROUTINE  DTIMECDRSEED, INPUT, OUTPUT) 
REAL *8  INPUT, OUTPUT, RV, DRSEED 
RV  =  GGUBFSC DRSEED) 

OUTPUT  =  LOG(RV)/( -INPUT) 

RETURN 

END 
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SUBROUTINE  DRGQTYCDRSEED, INPUT, OUTPUT) 
REAL *8  DRSEED, INPUT, OUTPUT, RV 
RV  =  GGUBFSC DRSEED) 

OUTPUT  =  AINTC  LOG( RV)/(LOGC INPUT) ) )+l 

RETURN 

END 


SUBROUTINE  CLASSC  DRSEED, INPUT, OUTPUT) 
REAL*8  DRSEED, INPUT, RV 
CHARACTER***  OUTPUT 
RV  =  GGUBFSC DRSEED) 

IFCRV  .IE.  INPUT)  THEN 
OUTPUT  =  'BAD' 

ELSE 

OUTPUT  =  'GOOD' 

ENDIF 

RETURN 

END 


SUBROUTINE  OPCOSTC DRSEED, INPUT1 , INPUT2, PI , P2, OUTPUT ) 
REAL*8  DRSEED, INPUT1 , INPUT2, PI , P2, OUTPUT, STOR, NSHIP 
OUTPUT  a  0 
JO  =  1/C1-P1) 

NSHIP  *  INPUT2*INPUT1 
OUTPUT  =  (P2*NSHIP)*J0 
RETURN 
END 


SUBROUTINE  RTIMEC DRSEED, INPUT1 , INPUT2, OUTPUT) 

C  EXPONENTIAL  DISTRIBUTION  -  BASE  MEAN  8  VARIANCE 
REAL*8  DRSEED, INPUT1 , INPUT2, OUTPUT, RV 
RV  =  GGUBFSC DRSEED) 

OUTPUT  =  L0GCRV)/C-INPUT1*INPUT2) 

RETURN 

END 


SUBROUTINE  RTIMEIC DRSEED, INPUT1 , INPUT2, OUTPUT) 

C  LOGNORMAL  DISTRIBUTION  -  BASE  MEAN  8  VARIANCE 

REAL*8  INPUT1,INPUT2, OUTPUT, STDV,AVG, PI, LOW, U1,U2,GG, 
8  DRSEED 

INTEGER  II, FINISH 
GG  =  LOGC2.0) 

LOW  =  999999 
FINISH  =  INTCINPUT2) 

STDEV  =  SQRTCGG) 

AVG  =  -CLOGCINPUT1))-  O.SXGG 
PI  =  3.141592654 
DO  662  II  =  1, FINISH 
U1  =  GGUBFSC DRSEED) 

U2  a  GGUBFSC DRSEED) 

RV  =  SQRTC-2*L0GCU1))*C0SC2*PI*U2)*STDEV+AVG 
RVT  a  EXPCRV) 

IFCRVT  .LT.  LOW)  THEN 
LOW  =  RVT 

ENDIF 

662  CONTINUE 

OUTPUT  *  LOW 

RETURN 

END 

SUBROUTINE  RTIME2CDRSEED, INPUT1 , INPUT2, OUTPUT) 

C  LOGNROMAL  DISTRIBUTION  -  BASE  MEAN  8  2XVARIANCE 

REALX8  INPUT1,INPUT2, OUTPUT, STDV, AVG, PI, LOW, U1,U2,GG, 
8  DRSEED- 
INTEGER  II, FINISH 
GG  =  LOGC3.0) 

LOW  »  999999 
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FINISH  =  INT ( INPUT2) 

STDEV  =  SQRTCGG) 

AVG  =  -( LOG( INPUT1) )-  0.5*GG 
PI  =  3.141592654 
DO  662  II  =  1, FINISH 
U1  =  GGUBFS(DRSEED) 

U2  =  GGUBFSC  DRSEED) 

R V  *  SORT C -2*L0G(U1) )*C0S(2*PI*U2)*$TDEV+AVG 
RVT  =  EXP(RV) 

IFCRVT  .LT.  LOW)  THEN 
LOW  =  RVT 

ENDIF 

662  CONTINUE 

OUTPUT  =  LOW 

RETURN 

END 

SUBROUTINE  RTIME3C DRSEED, INPUT1 ,  INPUT2,  OUTPUT) 

:  LOGNROMAL  DISTRIBUTION  -  BASE  MEAN  2  4*VARIANCE 

REAL *8  INPUT1, INPUT2, OUTPUT, STDV, AVG, PI, LOW, U1,U2,GG, 
2  DRSEED 

INTEGER  II, FINISH 
GG  =  LOGC5.0) 

LOW  =  999999 
FINISH  =  INT  < INPUT2) 

STDEV  =  SQRTCGG) 

AVG  =  -CLOG(INPUTl))-  0.5«GG 
PI  =  3.141592654 
DO  662  II  =  1, FINISH 
U1  =  GGUBFSC  DRSEED) 

U2  =  GGUBFSC  DRSEED) 

RV  =  SQRTC-2*L0G<U1))*C0S<2*PIXU2)XSTDEV+AVG 
RVT  =  EXPCRV) 

IFCRVT  .LT.  LOW)  THEN 
LOW  =  RVT 

ENDIF 

662  CONTINUE 

OUTPUT  =  LOW 

RETURN 

END 
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FILEi  TSTNUM  FORTRAN  A1 


THIS  PROGRAM  COMPUTES  THE  VALUES  OF  THE  THREE  REWARD  CRITERIA 
FOR  VARIOUS  VALUES  OF  L  FOR  COMPARISON  WITH  SIMULATION  RESULTS 

REAL  PB, LAMDA1 > LAMDA2. B ETA, MU, INC, ES, ED, FT, FTBAR, A, A BAR, 

8  L , TOEND, EC (100), EJ (100), EMC  100) , R1(100),R2(100),R3(100), 

8  ITIME(100),PG,B,BB,E,C1,C2,C3,C4,M1,M2,TEMP,J1,RATE 
INTEGER  MAXT , K 
MAXT  =  6 
INC  =  0.20 
PB  =  .20 
PG  =  1-PB 
LAMDA1  =1.10 
LAMDA2  =0.3 
ES=1/LAMDA1 
ED=1/LAMDA2 
BETA  =0.8 
MU  =  0.08 
B=BETA 
BB=1-B 

RATE  =  <LAMDA1*PB)/BB 
TOEND  =  MAXT/INC 
L  =  0 

DO  100  K  =  1 , TOEND 
L=L+INC 
£=EXP(-MUXL ) 

TEMP=l-(Bx£) 

FTBAR  =  (BBXEVTEMP 

FT=1-FTBAR 

ABAR=PBXFT 

A=1-ABAR 

Cl=  ES/ABAR 

C2=ED 

C3=CPB/ABAR)XC(BBX(BXMU))XL0GCTEMP/BB)-CLXFTBAR)) 

C4=(L*A)/ABAR 

Jl=( BB/FT)x( ( l/( BBXB8 ) )-(E/(TEMPXTEMP) ) ) 

M1=(C3+ED)XRATE+((LXAXRATE)/ABAR) 

M2=BBX(E/(TEMPXTEMP))X(1/ABAR)XPB 

ITIME(K)  =  L 

EC(K)  =  C1+C2+C3+C4 

EJ(K)  =  J 1 

EM(K)  =  M1+M2 

R1(K)  =  1/EC(K) 

R2(K)  =  EJ(K)/EC(K) 

R3(K)  =  (EJ(K)-EM(K))/EC(K) 

WRITE(2,132)  ITIME(K) ,R1(K),R2(K),R3(K) 

WRITE(3, 1321  ITIME(K) , EC(K), EJ(K) , EM(K) 

WRITE(4, *)  L 

WRITEC4, 133)  C1,C2,C3,C4 
WRITEC4. 134)  J1 
WRITE(4, 135)  Ml, M2 
FORMAT  ( IX, F6 . 2, 3F9 .3) 

FORMAT  (1X,4F12.5) 

FORMAT  ( IX, F12 . 5) 

FORMAT  (1X,2F12.5) 

WRITE(4,X) 

CONTINUE 

STOP  'PROGRAM  EXECUTION  COMPLETE' 

END 
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APPENDIX  B 


Recall  that  the  probability  distribution  for  the 
quantity  of  drugs  on  a  vessel  given  that  it  is  bad  is: 


bj  =  (1  -  3)(3J-‘;  P  =  0.8,  j  =  1,2,3, 


P{T>t| 


Ft  (t)  = 


Ft  (t)  =  S  (1-Fz  (t)  )->bj 

j-1 


1-PFz  (t) 


P { Tit }  =  1  - 


l-pST(t) 


1-Fz  (t) 
l-pfT(t) 


now  let 


1-PFz  (tP  ) 


which  can  be  written  as 


Fz  (tP )  * 


i-Pp 


Fz  (tp  )  -  1  -  Fz  <tp  )  =  pi1:!) 


Assuming  Fz  (t)  is  lognormal  with  mean  =  2.1792  and  standard 
deviation  =  0.8326  as  in  distribution  Case  A 


♦ 


P.(l-P) 

1-Pp 


therefore 


tp 


EXP 


U  +  (JZp  (  i  -  p  > 

l  -  0  P 


where  ZP  < i - o  >  is  taken  from  tables  of  the  Standard  Normal 

i  -  P  P 


distribution. 


Using  the  equation  for  tP ,  the  following  table  of  quantiles 
was  generated: 


p 

tp 

0.10 

1.65 

0.25 

2.46 

0.50 

3.95 

0.75 

6.78 

0.90 

11.99 

It  can  be  seen  that  these  quantile  values  agree  closely 
with  those  presented  in  the  sensitivity  analysis. 
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